ggcoverage User Guide
18 七月 2022
Abstract
The goal of ggcoverage is to simplify the process of visualizing genomic coverage. It contains functions to load genomic data with bam/bigwig/bedgraph file format, create genomic coverage plot, add various annotation to the coverage plot, including gene annotaion, transcription annotion, ideogram annotation and peak annotation.
Package
ggcoverage 0.5.0
1 Getting started
ggcoverage is an R package distributed as part of the Bioconductor project and CRAN. To install the package, start R and enter:
# install via Bioconductor
if (!requireNamespace("BiocManager", quietly=TRUE))
install.packages("BiocManager")
BiocManager::install("ggcoverage")
# install via CRAN
install.package("ggcoverage")
# install via Github
# install.package("remotes") #In case you have not installed it.
remotes::install_github("showteeth/ggcoverage")Once ggcoverage is installed, it can be loaded by the following command.
library("rtracklayer")
library("ggcoverage")2 Introduction
The goal of ggcoverage is simplify the process of visualizing genomic coverage. It contains three main parts:
- Load the data:
ggcoveragecan load bam, bigwig (.bw), bedgraph file from various NGS data, including RNA-seq, ChIP-seq, ATAC-seq, et al. - Create genomic coverage plot
- Add annotaions:
ggcoveragesupports four different annotaions:- gene annotaion: Visualize genomic coverage across whole gene
- transcription annotion: Visualize genomic coverage across different transcripts
- ideogram annotation: Visualize the region showing on whole chromosome
- peak annotation: Visualize genomic coverage and peak identified.
ggcoverage utilizes ggplot2 plotting system, so its usage is ggplot2-style!
3 RNA-seq data
3.1 Load the data
The RNA-seq data used here are from Transcription profiling by high throughput sequencing of HNRNPC knockdown and control HeLa cells, we select four sample to use as example: ERR127307_chr14, ERR127306_chr14, ERR127303_chr14, ERR127302_chr14, and all bam files are converted to bigwig file with deeptools.
Load metadata:
# load metadata
meta.file <- system.file("extdata", "RNA-seq", "meta_info.csv", package = "ggcoverage")
sample.meta = read.csv(meta.file)
sample.meta
#> SampleName Type Group
#> 1 ERR127302_chr14 KO_rep1 KO
#> 2 ERR127303_chr14 KO_rep2 KO
#> 3 ERR127306_chr14 WT_rep1 WT
#> 4 ERR127307_chr14 WT_rep2 WTLoad track files:
# track folder
track.folder = system.file("extdata", "RNA-seq", package = "ggcoverage")
# load bigwig file
track.df = LoadTrackFile(track.folder = track.folder, format = "bw",
meta.info = sample.meta)
# check data
head(track.df)
#> seqnames start end score Type Group
#> 1 chr14 21572751 21630650 0 KO_rep1 KO
#> 2 chr14 21630651 21630700 1 KO_rep1 KO
#> 3 chr14 21630701 21630800 4 KO_rep1 KO
#> 4 chr14 21630801 21657350 0 KO_rep1 KO
#> 5 chr14 21657351 21657450 1 KO_rep1 KO
#> 6 chr14 21657451 21663550 0 KO_rep1 KOPrepare mark region:
# create mark region
mark.region=data.frame(start=c(21678900,21732001,21737590),
end=c(21679900,21732400,21737650),
label=c("M1", "M2", "M3"))
# check data
mark.region
#> start end label
#> 1 21678900 21679900 M1
#> 2 21732001 21732400 M2
#> 3 21737590 21737650 M3Load GTF file:
gtf.file = system.file("extdata", "used_hg19.gtf", package = "ggcoverage")
gtf.gr = rtracklayer::import.gff(con = gtf.file, format = 'gtf')3.2 Basic coverage
basic.coverage = ggcoverage(data = track.df, color = "auto",
mark.region = mark.region, range.position = "out")
basic.coverage
You can also change Y axis style:
basic.coverage = ggcoverage(data = track.df, color = "auto",
mark.region = mark.region, range.position = "in")
basic.coverage
3.3 Add gene annotation
basic.coverage +
geom_gene(gtf.gr=gtf.gr)
3.4 Add transcript annotation
basic.coverage +
geom_transcript(gtf.gr=gtf.gr,label.vjust = 1.5)
3.5 Add ideogram
basic.coverage +
geom_gene(gtf.gr=gtf.gr) +
geom_ideogram(genome = "hg19",plot.space = 0)
#> [1] "hg19"
#> Loading ideogram...
#> Loading ranges...
#> Scale for 'x' is already present. Adding another scale for 'x', which will
#> replace the existing scale.
basic.coverage +
geom_transcript(gtf.gr=gtf.gr,label.vjust = 1.5) +
geom_ideogram(genome = "hg19",plot.space = 0)
#> [1] "hg19"
#> Loading ideogram...
#> Loading ranges...
#> Scale for 'x' is already present. Adding another scale for 'x', which will
#> replace the existing scale.
4 ChIP-seq data
The ChIP-seq data used here are from DiffBind, I select four sample to use as example: Chr18_MCF7_input, Chr18_MCF7_ER_1, Chr18_MCF7_ER_3, Chr18_MCF7_ER_2, and all bam files are converted to bigwig file with deeptools.
Create metadata:
# load metadata
sample.meta = data.frame(SampleName=c('Chr18_MCF7_ER_1','Chr18_MCF7_ER_2','Chr18_MCF7_ER_3','Chr18_MCF7_input'),
Type = c("MCF7_ER_1","MCF7_ER_2","MCF7_ER_3","MCF7_input"),
Group = c("IP", "IP", "IP", "Input"))
sample.meta
#> SampleName Type Group
#> 1 Chr18_MCF7_ER_1 MCF7_ER_1 IP
#> 2 Chr18_MCF7_ER_2 MCF7_ER_2 IP
#> 3 Chr18_MCF7_ER_3 MCF7_ER_3 IP
#> 4 Chr18_MCF7_input MCF7_input InputLoad track files:
# track folder
track.folder = system.file("extdata", "ChIP-seq", package = "ggcoverage")
# load bigwig file
track.df = LoadTrackFile(track.folder = track.folder, format = "bw",
meta.info = sample.meta)
# check data
head(track.df)
#> seqnames start end score Type Group
#> 1 chr18 76799701 76800000 439.316 MCF7_ER_1 IP
#> 2 chr18 76800001 76800300 658.974 MCF7_ER_1 IP
#> 3 chr18 76800301 76800600 219.658 MCF7_ER_1 IP
#> 4 chr18 76800601 76800900 658.974 MCF7_ER_1 IP
#> 5 chr18 76800901 76801200 0.000 MCF7_ER_1 IP
#> 6 chr18 76801201 76801500 219.658 MCF7_ER_1 IPPrepare mark region:
# create mark region
mark.region=data.frame(start=c(76822533),
end=c(76823743),
label=c("Promoter"))
# check data
mark.region
#> start end label
#> 1 76822533 76823743 Promoter4.1 Basic coverage
basic.coverage = ggcoverage(data = track.df, color = "auto", region = "chr18:76822285-76900000",
mark.region=mark.region, show.mark.label = FALSE)
basic.coverage
4.2 Add annotations
Add gene, ideogram and peak annotaions. To create peak annotaion, we first get consensus peaks with MSPC, you can also use DEbChIP’s GetConsensusPeak (MSPC’s wrapper) to do this.
# get consensus peak file
peak.file = system.file("extdata", "ChIP-seq", "consensus.peak", package = "ggcoverage")
basic.coverage +
geom_gene(gtf.gr=gtf.gr) +
geom_peak(bed.file = peak.file) +
geom_ideogram(genome = "hg19",plot.space = 0)
#> [1] "hg19"
#> Loading ideogram...
#> Loading ranges...
#> Scale for 'x' is already present. Adding another scale for 'x', which will
#> replace the existing scale.
5 Session info
sessionInfo()
#> R version 4.0.3 (2020-10-10)
#> Platform: x86_64-conda-linux-gnu (64-bit)
#> Running under: CentOS Linux 7 (Core)
#>
#> Matrix products: default
#> BLAS/LAPACK: /home/softwares/anaconda3/envs/r4.0/lib/libopenblasp-r0.3.12.so
#>
#> locale:
#> [1] LC_CTYPE=zh_CN.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=zh_CN.UTF-8 LC_COLLATE=zh_CN.UTF-8
#> [5] LC_MONETARY=zh_CN.UTF-8 LC_MESSAGES=zh_CN.UTF-8
#> [7] LC_PAPER=zh_CN.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=zh_CN.UTF-8 LC_IDENTIFICATION=C
#>
#> attached base packages:
#> [1] stats4 stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] ggcoverage_0.5.0 rtracklayer_1.50.0 GenomicRanges_1.42.0
#> [4] GenomeInfoDb_1.26.7 IRanges_2.24.1 S4Vectors_0.28.1
#> [7] BiocGenerics_0.42.0 knitr_1.37 BiocStyle_2.18.1
#>
#> loaded via a namespace (and not attached):
#> [1] colorspace_2.0-0 ellipsis_0.3.2
#> [3] biovizBase_1.38.0 htmlTable_2.4.0
#> [5] XVector_0.30.0 base64enc_0.1-3
#> [7] dichromat_2.0-0 rstudioapi_0.13
#> [9] farver_2.1.0 ggrepel_0.9.1
#> [11] bit64_4.0.5 AnnotationDbi_1.52.0
#> [13] fansi_0.4.2 xml2_1.3.2
#> [15] ggbio_1.38.0 splines_4.0.3
#> [17] ggh4x_0.2.1 cachem_1.0.4
#> [19] Formula_1.2-4 jsonlite_1.7.2
#> [21] Rsamtools_2.6.0 dbplyr_2.1.1
#> [23] cluster_2.1.1 png_0.1-7
#> [25] graph_1.68.0 BiocManager_1.30.16
#> [27] compiler_4.0.3 httr_1.4.2
#> [29] backports_1.2.1 assertthat_0.2.1
#> [31] Matrix_1.3-3 fastmap_1.1.0
#> [33] lazyeval_0.2.2 cli_3.3.0
#> [35] htmltools_0.5.2 prettyunits_1.1.1
#> [37] tools_4.0.3 gtable_0.3.0
#> [39] glue_1.6.2 GenomeInfoDbData_1.2.4
#> [41] reshape2_1.4.4 dplyr_1.0.5
#> [43] rappdirs_0.3.3 Rcpp_1.0.6
#> [45] Biobase_2.50.0 jquerylib_0.1.3
#> [47] vctrs_0.4.1 Biostrings_2.58.0
#> [49] xfun_0.30 stringr_1.4.0
#> [51] lifecycle_1.0.0 ensembldb_2.14.1
#> [53] XML_3.99-0.6 zlibbioc_1.36.0
#> [55] scales_1.1.1 BSgenome_1.58.0
#> [57] VariantAnnotation_1.36.0 ProtGenerics_1.22.0
#> [59] hms_1.0.0 MatrixGenerics_1.2.1
#> [61] RBGL_1.66.0 parallel_4.0.3
#> [63] SummarizedExperiment_1.20.0 AnnotationFilter_1.14.0
#> [65] RColorBrewer_1.1-2 curl_4.3
#> [67] yaml_2.2.1 memoise_2.0.0
#> [69] gridExtra_2.3 ggplot2_3.3.5
#> [71] sass_0.4.1 biomaRt_2.46.3
#> [73] rpart_4.1-15 reshape_0.8.8
#> [75] latticeExtra_0.6-29 stringi_1.5.3
#> [77] RSQLite_2.2.5 highr_0.8
#> [79] checkmate_2.0.0 GenomicFeatures_1.42.2
#> [81] BiocParallel_1.24.1 rlang_1.0.3
#> [83] pkgconfig_2.0.3 matrixStats_0.58.0
#> [85] bitops_1.0-6 evaluate_0.14
#> [87] lattice_0.20-45 purrr_0.3.4
#> [89] labeling_0.4.2 patchwork_1.0.0
#> [91] GenomicAlignments_1.26.0 htmlwidgets_1.5.3
#> [93] bit_4.0.4 tidyselect_1.1.0
#> [95] GGally_2.1.2 plyr_1.8.6
#> [97] magrittr_2.0.1 bookdown_0.26
#> [99] R6_2.5.0 generics_0.1.0
#> [101] Hmisc_4.6-0 DelayedArray_0.16.3
#> [103] DBI_1.1.1 pillar_1.5.1
#> [105] foreign_0.8-81 survival_3.2-10
#> [107] RCurl_1.98-1.3 nnet_7.3-16
#> [109] tibble_3.1.0 crayon_1.4.1
#> [111] utf8_1.2.1 OrganismDbi_1.32.0
#> [113] BiocFileCache_1.14.0 rmarkdown_2.14
#> [115] jpeg_0.1-8.1 progress_1.2.2
#> [117] grid_4.0.3 data.table_1.14.2
#> [119] blob_1.2.1 digest_0.6.27
#> [121] openssl_1.4.3 munsell_0.5.0
#> [123] bslib_0.3.1 askpass_1.1